A Lifting Relation from Macroscopic Variables 
to Mesoscopic Variables in Lattice Boltzmann 
Method: Derivation, Numerical Assessments 
and Coupling Computations Validation 

Hui Xu a ' b , Huibao Luan a , Yaling He a , Wenquan Tao a '* 

a Key Laboratory of Thermo-Fluid Science & Engineering, Xi'an Jiaotong 
University, Xi'an, Shaanxi, 710049, China. 

h Institut Jean le Rond d'Alembert, UMR CNRS 1190, Universite Pierre et Marie 
Curie - Paris 6, 4 Place Jussieu case 162 Tour 55-65, 75252 Paris Cedex 05, 

France 



* Corresponding author. Tel:+86-29-82669106. Fax:+86-29-82669106. 

Email addresses: xuhuixj@gmail.com or xu@lmm.jussieu.fr (Hui Xu), 

wqtao@mail.xjtu.edu.cn (Wenquan Tao ). 



Preprint submitted to Elsevier 20 January 2013 



Abstract 



In this paper, analytic relations between the macroscopic variables and the 
mesoscopic variables are derived for lattice Boltzmann methods (LBM). The 
analytic relations are achieved by two different methods for the exchange from 
velocity fields of finite-type methods to the single particle distribution func- 
tions of LBM. The numerical errors of reconstructing the single particle distri- 
bution functions and the non-equilibrium distribution function by macroscopic 
fields are investigated. Results show that their accuracy is better than the ex- 
isting ones. The proposed reconstruction operator has been used to implement 
the coupling computations of LBM and macro-numerical methods of FVM. 
The lid-driven cavity flow is chosen to carry out the coupling computations 
based on the numerical strategies of domain decomposition methods (DDM). 
The numerical results show that the proposed lifting relations are accurate 
and robust. 

Key words: LBM, Navier-Stokes equations, non-equilibrium distribution func- 
tions, multi-scale perturbation expansion, coupling computation, FVM 
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1 Introduction 



In the past decades, LBM has been widely used to simulate fluid flow problems 
[Tf2] , including complex turbulent fluid flows [3|H] and multiscale modeling [5|6] . 
This method is based on the Boltzmann kinetic equation which is used to de- 
scribe a number of interacting populations of particles. As described in [7], 
"The LBE could potentially play a twofold function-as a telescope for the 
atomistic scale and a microscope for the macroscopic scale". In [S] dense flu- 
ids flow past and through a carbon nano tube (CNT) was studied by a hybrid 
model coupling LBM and MDS. The authors pointed out that replacing the 
finite volume solver by a LBM aims to take advantage of the mesoscopic 
modeling inherent in LB simulations. Thus LBM is a mesoscopic method in 
nature is a widely-accepted understanding in the literature. The macroscopic 
parameters such as fluid density, velocity and pressure can be obtained via 
some averages of the mesoscopic variable which conform the basic conservation 
laws of mass and momentum [2]. In practical applications of LBM to simulate 
a macroscopic problem, a crucial problem is confronted, that is, a reasonable 
initial meso-field must be specified to start the evolution process. The first ini- 
tializing method was proposed in [9] in 1993. Recently, several methods have 
been proposed to improve the accuracy of numerical results and reduce the 
initial layers (oscillation layers) [TUfTY] . Such oscillations have a numerical ori- 
gin and are due to the artificial compressibility of LBM. Here, " initial layer " 
refers to such a computational stage within which the macroscopic parameters 
are oscillating. When the initial data is not well-prepared, there is an initial 
layer during which the solution adapts itself to match the profile dictated by 
the environment. For the LBM, the existence of the initial layers is a common 
phenomenon [10]. In this paper, we will derive the lifting relations between 
the macroscopic variables and the mesoscopic variables in LBM by two ways. 
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According to the authors' knowledge, the proposed lifting relations in this 
paper are different from those in the existing literature PllUHllfl2fl3H14fl5j . 
The proposed relations will offer us some new views about the reconstruction 
of nonequilibrium distribution functions in LBM. 



Challenging multiscale phenomena or processes are widely existed in material 
science, chemical engineering process, energy and power engineering, and other 
engineering fields. Generally speaking, for a multiscale problem, we often must 
use different methods to numerically model the processes at different geomet- 
ric sub- regions and exchange solution information at interface ^6|17IH8|IT9] . 
Such coupling computations are widely adopted in the present-day multiscale 
simulation. As indicated above LBM is a kind of mesoscopic methods, which 
is a candidate to implement the meso-macro or micro-meso coupling compu- 
tations in engineering applications [TJ. So, the proposed method not only can 
be used to obtain a better initial field for LBM, but also can be adopted in 
the multi-scale computation. For example in [7J the possibility of coupling 
LBM with molecular dynamics simulation (MDS) was investigated and found 
that with proper time and geometric scales the two numerical methods can 
be coupled. And in [8] such coupling simulation was conducted. In the ex- 
isting literatures the coupling of finite difference method (FDM, which is a 
macrosopic method) with LBM was adopted in [T9"|2Uf21| . but the proposed 
coupling method is similar to a multigrid method and a simple regulariza- 
tion formula is used in their computations. The regularization formula in [19] 
only considers the first-order approximation of the single particle distribu- 
tion function and the coupling formula in [20] is only used to deal with the 
one-dimensional reaction-diffusion system. In [8] the coupling between LBM 
and MDS was implemented by exchange of velocity and velocity gradient at 
the interface region. In this paper, the proposed meso-macro (or micro-meso) 
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coupling is expected to be used for domain decomposition methods, in which 
LBM and macro-type numerical method (or micro-type numerical method and 
LBM) are adopted in different sub-domain and information is exchanged at the 
interface. We believe that our proposed relation is more useful method for en- 
gineering multiscale computations. In addition, the proposed coupling method 
can also be used to carry out the multigrid computations and equation-free 
multiscale (EFM) computations [22]. It is well-known that LBM is very power- 
full for the parallel computing on a low cost [23f2l] . So, the proposed relation 
can be used in the parallel simulations for multiscale simulations of complex 
fluid flows based on the refinement strategies. 

To the authors' understanding the glossary "lifting relation" means that macro- 
scopic variables in a lower degree-of-freedom (DoF) system are upscaled to 
meso/microscopic variables in a higher DoF system. Generally, it is difficult 
to establish the one-to-one map from a lower DoF system to a higher DoF 
system, although the lower DoF system can be seemed to be an approximate 
or approaching form of a higher DoF system in some referred scales. This sit- 
uation happens when numerical results of different scales are coupled at the 
same location. For example when MDS and continuum method are coupled, 
reference [25] indicated that it is straightforward to obtain the continuum 
quantities (such as velocity, pressure) from the particle description by averag- 
ing over the local region and over time, but the reverse problem, generating 
meso/microscopic particle configuration from known macroscopic quantities is 
non-trivial and must necessarily be non-unique. The glossary "lifting relation" 
in the title of this paper is proposed based on the concept of the DoF of the 
governing equations. 

In this paper, we will give two methods to establish the relations between 
variables of the Navier-Stokes equations and variables of LBM. Numerical 
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tests demonstrate that the proposed methods of computing non-equilibrium 
distribution functions are effective and accurate. 

The rest of the paper is organized as follow. In section |2j the details of multi- 
scale derivation of non-equilibrium distribution functions is given. In section |3j 
the non-equilibrium distribution functions are obtained by Boltzmann-BGK 
equations. In section |4j the performances of the proposed relations to recon- 
struct non-equilibrium distribution functions are demonstrated by numerical 
tests. Finally, some conclusions are given. 

2 Lattice Boltzmann hydrodynamics and multiscale approach 

In this section, we will review LBM and the corresponding macroscopic equa- 
tion. Based on this review, we will derive a relation for lifting macroscopic 
variables to microscopic variables by multiscale approach. 

2.1 Lattice Boltzmann hydrodynamics 

We now introduce the lattice Boltzmann-BGK model as a solver for the 
weakly-compressible Navier-Stokes equations. LBM is built up from the lattice 
gas cellular automata models [2]. The numerical scheme of LBM is established 
based on a finite discrete- velocity model of the Boltzmann-BGK equation and 
can be expressed as follows 

/ i (x + (Jtc i ,t + <Jt)-/(x J t) = n i , (i) 

where fi represents the single-particle distribution function along the direction 
Cj ( % — 0, . . . , n), Cj is the element of the discrete velocity set V = {c , . . . , c n }. 
Qi denotes the collision operator which is non-dimensional. The macroscopic 
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variables, the density p and the velocity u, are defined locally by the distribu- 
tion functions as follows 



pM) = E/*M) = E/rM), (2) 

j=0 i=Q 

u(x,f) = 1 E */*(*,*) = 1 E cji° q) (x,t). (3) 

For the standard LBM, the collision operator is defined by the so-called BGK 
collision 



BGK 



— m^t)-ft q \x,t)]. (4) 



n 

Tlbr 

For the convenience of comparison, from here, we use the similar notations in 
[26] . The local equilibrium distribution is defined by 

ft q \x,t) = f l Lieq) (x,t) + f^\x,t), (5) 

where (x, t) and (x, t) denote the linear part and the quadratic part 
of the equilibrium distribution, respectively. The linear part is given by 



^ (eq) (x,t)=^p(l + ^c,-u(x,t)), (6) 

and the quadratic part is expressed by 

/? (eq) (x,t) = ^(MWM)) : Ei, (7) 

where c s is the lattice sound speed of the model, Ui denotes the weight and 
Ej is a second-order tensor defined by 

Ej a/ g = CiaCip — (? s 5 a p. (8) 

The tensor product definition between two first order tensors a and b is given 
as follows 

(ab) Q/3 = a a bp, (9) 



1 
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and the corresponding second-order tensor :-product between A and B is given 
by 

d 

A : B = ]T A «/3 B */3, (10) 

a,P=l 

where d denotes the spatial dimension. 

In this paper, we mainly focus on the standard LBM. By the Chapman-Enskog 
expansion, under the small Ma number restriction (Ma < 0.2), we can recover 
the Navier-Stokes equations as follows 

d tP + d a (pu a ) + 0(5t 2 ) = 0, (11) 

d t (pu a ) + dp(pu a up) = -d a p + ud p (p(d a up + d^u a )) + 0(8t 2 )+0(5tu 3 ), (12) 
where p is defined by 

P = c 2 s p- 

It is clear that the recovered Navier-Stokes equations are weakly compressible 
[211271 28] . So, the density is coupled with the pressure field in LBM. In Eq. 



(12), the second term of R.H.S can be rewritten as 

vdp(p(d a up + dpu a )) = vp{dpdpUa) + v(dpp)(d a up + d p u a ) + vpd a dpu p . (13) 
And the corresponding third-order term 0(Stu 3 ) is given by 

0(5tu 3 ) = —adpd^pUaUpUy). (14) 
The fluid viscosity v is defined by 

v = c 2 s (ri hm - -)6t, (15) 



and a is given by 



a = - 9 . (16) 



In Eq. (13), the third term of R.H.S will vanish for a divergence-free field. 



But the second term will not vanish, if the density p is nonhomogeneous in 
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the spatial domain. The Navier-Stokes equations are recovered by LBM under 
the low Mach condition. Physically, LBM is a weakly compressible model for 
solving Navier-Stokes equations. 



At this point, we describe two situations where the lifting relation is useful. The 
first situation is using the lifting relation to get a good initial field of the density 
distribution function from specified velocity and pressure fields. As indicated 
above the recovered Navier-Stokes equations are weakly compressible, hence 
pressure field is coupled with the density field by the equation of state (p = 
cj.p). In engineering computations, the weakly-compressible flow is often used 
as an approximation of the incompressible flow. For the lifting function, the 
consideration should be made from the weakly compressible side. The non- 
homogeneous character of the initial density is very significant for an initial 
routine of LBM in the proposed lift relation. This significance can be observed 
from the follow-up derivations. For the initial processes, if the initial pressure 
field is given, the lifting relation can be used to obtain the initial distribution 
functions consistent with the recovered Navier-Stokes equations. In another 
development when we couple LBM with other macroscopic solver of Navier- 
Stokes equations, we need to pass the macroscopic variables (pressure and 
velocity fields) to an approximate single particle distribution functions or the 
non-equilibrium distribution functions. At this time, a macroscopic equation 
relating to the given velocity and pressure to the particle distribution function 
of LBM become very useful. The major goal of the present paper is to derive 
such a lift relation, or a reconstruction operator as depicted in [9]. 

For the convenience of deriving such an equation, some changes are made for 



the form of Eq. (12). We first rewrite Eq. (12) as 



d t (pu a ) + dp(pu a up) = p{d t u a + updpUa) + u a (d t p + dp{pup)). ( 17 ) 
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If the initial velocity field is divergence- free, we have 



d t (pu a ) + dp(pu a u p ) = p(d t u a + dfj(u a up)) + u a (d t p + dp(pup)). (U 



The neglecting of the term pu a dpup is a widely accepted approximation. Ac- 



cording to Eq. (11), we have 



d t (pu a ) + dp{pu a u p ) = p(d t u a + dp{u a u p )). 



(19) 



Now, combining Eq. (12), Eq. (13) and Eq. (17), we gain 



d t u a + updpUa = -&f + v(dpdpu a + d a d p up) + v^f-(d a u p + dpu a ). ( 20 ) 



2.2 Derivation of Non-equilibrium Distribution Function by Multi-scale Ap- 
proach 



The coupled macro- micro/mesoscale simulation is a rapidly developing area 
of research that deals with processes covering several order of geometries. For 
such numerical approach, one needs to construct an initial condition u(x, 0) for 
the meso/microscopic simulator, which is corresponding to the initial macro- 
scopic variable U(x,0). Here, u(x,0) represents the meso/microscopic state 
variables and U(x, 0) stands for macroscopic state variables. As indicated 
above this procedure is called lifting [22] or reconstruction [29] step. The 
lifting (reconstruction) operator p is defined by 

u (x,0) = p{U(x,0)). (21) 

The lifting procedure leads to a one-to-many mapping. After the initialization 
of the meso/microscopic variables by the reconstruction operator p, they will 
be evolved by the meso/microscopic simulator. In this paper, LBM is adopted 
as the mesoscopic simulator. As indicated in [T51I2T)] the macroscopic state 
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variables are easy to be achieved. To transfer the micro/meso-scale parameters 
into macro parameters we need some restriction [22] or compression [2H] 
operators. Conceptually, this operator Ai is defined by 

U(x,t) = M{u{x,t)). (22) 

For LBM, the operator Ai is implemented by Eq. ^ and Eq. ([3]). Our at- 
tention will put on the development of the reconstruction operator /i by the 
multi-scale analysis. As discussed above the reconstruction operator in multi- 
scale computation is corresponding to the lifting relation in an initial problem. 
In the following we will derive the operator from the initial problem aspect. 

To obtain an appropriate initial field, we turn to a simple multiscale pertur- 
bation expansion. We separate the time scale into two different time scales, 
t\ = et (diffusive time-scale) and t 2 = e 2 t (convective time-scale). The time 
derivative dt is expanded using a small parameter e, which normally is pro- 
portional to the small Knudsen number (Kn < 0.1) [27J, 

d t = ed tl +e 2 d t2 + 0(e 3 ). (23) 

Similarly, introducing space scale x\ = ex, the corresponding spatial derivative 
is not expanded beyond the first-order term [27] 

d a = ed la + Q(e 2 ). (24) 

The single-particle distribution function is expanded as follows [27] 

ffat) = // 0) (x,t) + e /, (1) (x,f) + e'/fW) + • . . . (25) 

By the Taylor expansion, from Eq. (fi]), we get 

5t(d t + cMf^t) + 5t 2 (d t + c %a d a ) 2 f % {^t) + 0(5t 3 ) = SV (26) 
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Combining Eq.(23)-Eq.(25) with Eq.(26), we obtain 



(eq)/ 



(27) 



and 



efPfot) + e 2 tf>(x,t) = -r lhm [(ed tl + e 2 d t , 2 + ec ia d la )5t+ 

{\e 2 d 2 tl c ia d la + e 2 d h c ia d la + \e 2 c ia c ij3 d 1(X dxp)8t 2 ] 

(^ (0) (x,t) + 6/i 1) (x,t)) + O(5t 3 ). 
For first order of e, we get 



,2.(2) 



if } (x,i) = -n bm St{d tl + c ia d la (x,t))fP> + 0(6t d ) 



(cq) 



(2* 



(29) 



According to Eq. (|2])-Eq. ([3]), we have following equations in the first-order 
scale of e [30] 

d tlP + d la (pu a ) + O(5t 2 ) = 0, (30) 
d tl {pu a ) + d X p(pu a u p + c 2 s p5 aP ) + 0(5t 2 ) = 0, (31) 



Then, Eq. (31) can be rewritten as 



pd tl (u a ) + pUfid w (pu a + c 2 s p5 al3 ) + 0(5t 2 ) = 0. 



(32) 



By matching small scales, from Eq. (28), we can get up to the second order 
equations of the small parameter e: 



/f = -Ti hm 5td t jr - 5t 2 (r lhm - t){d tl + c^pffP + 0(5t 3 ). (33) 



e(0) 



1 



(34) 
(35) 



Then, we can get [30] 

d t2 p + O(5t 2 ) = 0, 
d t2 (pu a ) = vdip(p(di a up + dipu a )) + 0(5t 2 + 8tu 3 ) 



Furthermore, from Eq. (29), we have 



/ 4 (1) (M) = -r lhm 5t(d tl + c ia d la )(fr q, (x,t) + tf» w (x,f)) + O(^). (36) 



cL(eq) 



Q(eq) , 
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In the derivation of Eq. (36), we introduce the following formulas according 
to the chain rule of derivatives |3T] 



d t jr\x,t) = dJ e *\x,t)d tl p + d u J e(! \x,t)d tl up, (37) 

«9 lQ /f q) (x, t) = d p ft q) (x, t)d laP + d U0 f^\x, t)d la u p . (38) 
Now, the equilibrium function can be differentiated by the macroscopic vari- 
ables as follows [21] 

d P /i eq W) = ^M), (39) 
<M eq) (x, t) = d u J^\x, t) + cV« (eq) (x, t). (40) 



According to Eq. ^ and Eq. ([7]), we have 



du p fl [eci) = Uip^c il3 , (41) 



1 



So, we have 

d Ul3 ft q) = Uip[\{cip - up) + \c ia Cipu a ]. (43) 

C s C s 



Come here we can have following corollaries. 



Corollary 1 From Eq. (36), for the first-order approximation of e, there ex- 
ists a lifting relation from the macroscopic variables to the microscopic variable 

fi = -n hm St{(c ia - U a )di a pd p f- Q + (C ia - tO^a«00«p/i - 

pd p ft q) d a u a - - p d la pd Ua ft q) } 
= -Tibm5t{(c ia - u a )- p d la pfl e(i) + (c ia - u a ) ( 44 ) 

di a UpUip[~s{Cip - Up) + ^C i/ gCi 7 M 7 ]- 

ft q) d la u a - -dipp uJip[js{cip - up) + 4jc i/3 Ci 7 M 7 ]}. 
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Corollary 2 From Eq. (33), for the second-order scale of e, we have the fol- 
lowing approximation 

fP « -r lhm 5td t J C(1 \ (45) 
where the second-order derivative of is ignored. 



(2) 

Hence, we can easily establish an approximation for f- by the method anal- 
ogous to the approximation of as follows 



dtjt q) (*, t) = d p fP>(x, t)d t2 p + d up fr>(*, t)d t2 u p . 



e(eq) 



{eq), 



(46) 



By Eq. (34), we have 



d t2 fi eq) (^,t) = d u jP'(x,t)d t2 up = -d^fP 1 {x,t)d t2 (pup). 



r(eq)/ 



c(eq)/ 



From Eq. (35) and Eq. (43), it is easy to obtain 



(47) 



dt 2 fi Cq) = vui[^(ci/3 - up) + ^c i pc il u 1 }d la (p(d 1 pu a + d la up)). 

C s C s 



So, we have 



(4£ 



e 2 // 2) ~ -TdtuUi[^(Cip - Up) + ^ i C i pCryU- y }d a (p(dpU a + 3 a Up)). (49) 

By a simple derivation, we have 



1 



d a (p(dpu a + 9 a «^)) = d a p(dpu a + d a up) + p(dpd a u a + d 2 a up). (50) 



From Eqs. (49)~(50), we have 



eVi ~ -Tlbm^^i[^(Ci/3 - Up) + jiCipC^U^dapidpUa + O a Up) + 

p(dpd a u a + d 2 a up)) 



(51) 



Therefore, we get the following approximation of the non-equilibrium distri- 



bution function from Eq. ( 25 ) 



(52) 
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that is, 



fj neq \x,t) ft! -ri bm 5t{M T ,ia^ a p/f q) + U T ,iad a Ul3Uip[jlUT,i/3 + ^CiflC^U,,} 

+u}ip[jsUT,i/3 + jiC il3 Ci 7 u 7 ](v~d a p(dpu a + d a u p ) + v{dpd a u a + dfcup))}, 

(53) 

where Ur,i a = Cia — u a (ut = c« — u, peculiar velocity). Since the velocity field 
is divergence-free, we have 

^(neq.dfv) _ - Tlbm (ft i<9 a p// eq) + U T ,i a d a UpUJip[^UT,ip + ^CipC^Uj}- 

+u)ip[jzU T ,ii3 + ^CipCi 7 u y )(f^d a p(d^u a + d a u p ) + vd 2 a up)}. 

(54) 

Here, we also introduce an approximation of d Up f- by ignoring the higher- 
order terms of u 2 as adopted in [3T] 



dupft q) = Uip[\u T ,ij3 + ft ^£/ f (eq) . (55) 



Now, we have 



y(neq-dfv) _ -Ti hm 8t^ f-^Ur ,i(3 (u T ,i a d a U/3 + V-d a p{dpU a + d^) + vd\up). 

(56) 

Rewriting the above formula, we obtain 

^ neq _dfv) _ -T lhm 5tf^ od -^u T;i p{u Tiia d a up + vd 2 a up + v±d a pS af} ), ( 57 ) 
where = dpUa + d a up. 

In all, we can get an approximation of the single-particle distribution function 
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for divergence-free velocity fields as follows 

ft ~ ft^i 1 ~ ^Ti hm 6tu T> ifi(u T ,iadaUl3 + vd 2 a Up + U^d a pS a ^)}, ( 58 ) 

By a similar deviation, we can get an approximation of the single-particle 
distribution function for weak-compressible velocity fields as follows: 



fi ~ /i - ^nbm6t[u T ,ip(UT,iad a Up + v{p 2 a Ufi + dpd a U a ) 

+v l d a pS a p) - c 2 s d a u a )}, 



(59) 



Now we compare our results with that published in literatures. 



1. T. Imamura et al^T\ obtained the following formula 

/^/i eq) + e/i 1) = /i cq) [l-r lbm( 5t( 



<«0 , _ # («0n ... ^ u T, ia u T , _ _ (6Q) 



They only used to approximate the single-partial distribution functions. It 

(2) 

is well-known that in order to recover the correct Navier-Stokes equations, f} 



is needed. From this point of view, the approaching form (59) of the distribu- 



tion functions are more accurate than (60). If the divergence- free velocity field 



is considered, Eq. (58) is also superior to Eq. (60) because Eq. (58) contains 



(2) 

the information of fl which is related with molecule viscosity and density 
gradient. As for the lifting relation it is certainly essential to involve molecule 
viscosity and density gradient [T9p0] . 

2. Skodors [H] gave the following formula (ignoring the term of 0(Ma 2 )) 

1 



^neq),S = 



r CjCj : V(pu) - V ■ (pu) 



(61) 



Guo and Zhao [32] further simplified Eq.(61) and obtained the following rela- 
tion 



f 



(neq),G <- PO ^ 

= -TlbnA^i^CjCj : VU. 
c 



(62) 
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It is very clear that Eqs.(58) and (59) are totally different from Eqs.(61) and 



(62), respectively. The co-existence of Eqs. (58)~(59) and Eqs.(61 )~(62) as 
the lift correlation for the same situation may be regarded as the witness that 
the transformation from one-to-many must necessarily be not unique [25] . 



Some comparisons will be performed in Sec. 4 between Eqs. (58)~(59) and 



Eqs. (61 )^62) for schemes of D2Q9 and D2Q17. It turns out that the accuracy 



of Eqs. (58)~(59) derived in this paper is better than that of Eqs. (61 )~(62 ). 



(2) 



The derivation procedures of Eqs. (58)~(59) kept the information of the /. 
and other more details which are important to reduce the reconstruction rel- 
ative errors. 



3 Derivation of Non-equilibrium distributions via Boltzmann-BGK 
equations 



The Boltzmann equation [33] describes the statistical distribution of particles 
in a fluid. It is one of the most important equations of non-equilibrium statisti- 
cal mechanics, which deals with systems far from thermodynamic equilibrium 
j. The Boltzmann equation is described by 

d /( *; V,t) + v • V,/(x, v, t) + ±F(x) V,/(x, v, t) = fi(/(x, v, t)). (63) 



The Boltzmann equation (63) is an equation for the time t evolution of the 



distribution (properly a density) function f(x, v, t) in one-particle phase space, 
where x = (xi,X2, ■■■ ,Xd) £ R d and v = (i>i,i>2, ••• ,Vd) € -R d (d denotes 
the spatial dimension) are position and velocity, respectively. The equilibrium 
distribution function f eq (x,v,t) can be determined by 



rn 



27r«T(x,t)_ 



cxp 



ni 



2«T(x,t) 



u{x,t)Y 



(64) 
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Here, the quantities T(x,t), n(x,t) and u(x, t) represent the local tempera- 
ture, the local particle-number distribution density and the local velocity [2"||3"1] , 
repectively. ut = v — u(x, t) is the so called thermal velocity, m denotes 
the single-particle mass which is set to be unity for convenience. In order 
to simplify the complex collisional term, the following conserved relaxation 
time approximation is used to describe the collision term through only one 
characteristic frequency [31] 

d /( g t V,t) + v ■ V/(x, v, t) = -i(/(x, v, t) - /(«0(x, v, t)), (65) 

where the external force term is not considered and V denotes V a . r represents 
the relaxation time. 



In order to solve Eq.(65), the velocity space is discretized [2] and we gain 



+ Cl • V/,(x,f) = ~{Mn,t) - // eq) (x,*)), (66) 

where W{ denotes the integral weight factor, /j(x, t) = Wif(x, c«, t) and f- eq \x, t) = 
Wif^(x, Cj, t). Furthermore, along the characteristic line, the time-discretization 



form of Eq.(66) can be expressed as [2,36 

1 



fi(yL + Ci6t,t + 6t) = fi(x,t) - —(f l (x,t)-ft q \x,t)), i = 0,1,..., iV. (67) 

Tlbm 

where fi is the probability distribution function (PDF) along the ith direction, 
// eq ' ) is its corresponding equilibrium PDF, 8t is the time step, Cj is the particle 
velocity in the ith direction, and N is the number of the discrete particle 
velocities. Note: rib m = r/8t which is a dimensionless relaxation time. The 
local macro quantities are defined by Eqs. @ and Q. 

At the low fluid flow velocity (or low Mach number), an approximate form of 
the equilibrium distribution function fj^ is described by the discrete equilib- 
rium distribution, Eqs. ([5])~([7]). 
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Now, we consider the conserved relaxation time approximation of the Boltz- 



mann equation (65). The right hand side of Eq. (65) represents the difference 



between the distribution function and a local Maxwell distribution. This dif- 
ference is termed non-equilibrium distribution defined by 



/ (ncq) (x,v,t) = /(x,v,t) - /^(x,v,t). 



(eq)f 



(6£ 



Then, Eq. (65) can be rewritten as follows 



{wt +v ' v ) /(neq) ( x ' v '*) + (I + v ' v ) / (eq W,0 = -^/ (ncq) (x,v,t). 

(69) 

In the hydrodynamic region [33], the first term on the left-hand side of Eq. 



(69) can be neglected compared with the right-hand side [34] . Then, we obtain 



' d 
dt 



•V) f^\x 1 v 1 t) = -^f^\x,v,t). 



(70) 



In terms of the Maxwell equilibrium distribution and assuming a uniform 
temperature of the system, we can obtain 



/ (eq) (x,v,t) / d 



n(x,t) 



(| + v • V) n(x, t) - (x, v,t) (| + v-V)| 

= _I / ( M q)( XjVjf ) j 



2 
T 

2kT 



(71) 



where T = T(x, t) = constant. In Eq. (71 ), the left-hand term can be rewritten 
as follows 



(|+v-V)n(x,t) = (| + u(x,t).V)n(x,«) 



n(x,t) 



(72) 



+ n(xj) U T- Vn( X ,t) 



In order to satisfy the mass conservation condition of the fluid flow system, 



the first term of the right-hand side in Eq. (72) should be equal to zero. Hence 
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we have the following equation 



/<«0(x, v, t) (| + v • V) § - ^^f^ux • Vn(x, t) - n(x, t) V ■ u(x, f)) = 

i/^(x,v,t), 

(73) 

where c s = a/ kT. The term up is the thermal fluctuation energy, thus the 
non-equilibrium is determined by the material derivative of this thermal fluc- 
tuation energy. The quantity ( J| + v ■ v) u|> can be determined by the dynam- 
ical equation corresponding to the micro dynamical system. Here, we rewrite 
I (§i + v ' ^) u t as follows 



Hit + V ' V ) U ^ = ~ UT ' (It + U ^ X ' ^ ' V ) U(X ' ^)- u t- (ut-V)u(x, t). (74) 



Generally, the governing equation of the macroscopic physical quantity is rep- 
resented by 



^u(x, t)=(^ + u(x, t) ■ V j u(x, t) = F(x, u(x, t),t). (75) 



Normally, the macroscopic physical quantity u(x, t) in the governing equation 
is known. So, F(x, u(x, t),t) can be determined easily. For fluid flow problems, 
taking u(x, t) as fluid velocity, then F(x, u(x, £),£) can be estimated by fluid 



acceleration. The term up • (up • V)u in Eq. (74) can be determined by u(x, t) 
and the spatial derivatives of u(x, t). 

The lattice Boltzmann model is a special discrete form of the BGK lattice 
Bolzmann equation with respect to temporal and spatial variables. For LBM 
the equilibrium distribution, Eq. (|5]), is a polynomial-truncated approximation 



of the Maxwell distribution up to 0(|u| 3 ), so Eq. (73) can be applied to LBM 
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directly as follows 



// eq) (x,t) (o + Ci • V) % - ^gy%i >T " Vn(x,t) - n(x,t)V ■ u(x,t)) 



(76) 



where u^t = Cj — u(x, £). Now, the non-equilibrium distribution function can 
be denoted by 



/? eq) M) 



r/| eq) (x,t) 



u ljT • (F(x, u(x, t),t) + (ui jT • V)u + ^ Vn(x, t)J - c^V • u(x, t) 



(77) 



The derivation of Eq.(77) is completed based on the rigorous inherent physical 



consistency in the hydrodynamic region and and the derivation is independent 
on the spatial dimension. Meanwhile, the Maxwell equilibrium distribution is 
regarded as the tool to implement the analysis. 

It is worth pointing out that for DnQb LBM, F(x, u(x, t), t) can easily be 
determined from the recovered Naiver-Stokes equations, so the obtained non- 



equilibrium distribution function formulas (77) and (59) are identical. Thus 



by using different derivation method we come to the same conclusion. 



In addition, according to Eqs. (58), (59) and (77), it can be seen that the 



non-equilibrium distribution functions have the following form 



/i neq) = /i eq) A,(p,u), 



where \i(p, u) is a perturbative parameter with respect to p and u. The pa- 



rameter Aj(p, u) in Eq. (78) needs to satisfy the following constraints 



E/f q) A l (p,u) = 0, $>/r j A,(p,u) = 

i c,eV 



(79) 
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4 Numerical Tests 



In this section, the non-equilibrium distribution function will be validated by 
numerical tests. The numerical tests focus on validating the precision of the 
reconstruction operator and the correctness of the coupling computations. It's 
worth noting that the word "multiscale simulation" used in this paper is re- 
ferred to the coupling between numerical methods of microscale (molecular 
dynamics simulation), mesoscale (LBM) and macroscale(say, FVM) adopted 
in neighboring computational regions. And for such coupling the major con- 
cern is the transformation of solutions from macro(or meso)scales to meso(or 
micro)scales at the interface. The focus of the following presentation is to val- 
idate the correctness of the proposed operators. Because of space limitation 
the effect of the grid fineness on the numerical solution will not be conducted. 
Reference [33] can be referred. The effect of the mesh size on the accuracy of 



the reconstruction operator will be presented in Sec. 4.2 



4-1 Examination of the precision of the reconstruction operator 



In order to validate Formula ([77]), the D2Q9 [36J and D2Q17 [37] LBM are 
adopted to simulate 2D fluid flows. At low Mach number (Ma = u(x, t)/c s <C 



1), the R.H.S of Eq.(75) is equal to the R.H.S of Eq. (20) 



F a (x,u(x, t),t) = -— + v{dpdpu a + d a dpUp) + v^-{d a up + d^u a )(80) 



where 



V = C^(Tib m - -)St. 



51) 



The details of the macroscopic dynamic equation corresponding to D2Q17 



LBM are omitted (see [37]). Now, the non-equilibrium distribution in Eq.(77) 
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can be determined directly by the right-hand side of Eq.(20). For any given 
initial velocity and density fields, each term in the right-hand side of Eq.(20) 
can be calculated. In order to validate the precision of the proposed method, 
the following two basic quantities are defined 

f^t) = ft q \x,t) + fi^\x,t), (82) 

/f nCq) (x,t) = / J (x,t)-# q) (x,t) (83) 
where f} neq ^ (x, t) is called reconstructed non-equilibrium distribution function 



and is calculated by Eq.(77) and /j(x, t) is the reconstructed single-particle 
distribution function. f- r (x, t) and /j(x, t) denote the real non-equilibrium 
distribution function and the real single-particle distribution function, respec- 
tively. Here, we give two kinds of relative error definitions: single particle 
distribution function reconstruction error, single particle non-equilibrium dis- 
tribution function reconstruction error 



1 |/ t ( X ,t)-/,(x,t)|2 

\Numx(n + l) 2 ? 2 t ' ffat)* ' 1 J 



j^y(neq) ytneq)x 



i EE l/7 ncq W)-/f cq) M)l 2 (g5) 



" \ Num x (n + 1) V Y if ° q) (x, i) 2 

where Num denotes the number of lattice nodes. 

In order to demonstrate the proposed method, a freely-decaying 2D turbulence 
problem will be simulated by the proposed method. This turbulence problem 
often makes the local discrete single-particle distribution functions to be far 
from the local discrete equilibrium distribution functions, which yields a rich 
velocity structure. The freely-decaying 2D turbulence is implemented in a 
periodic box Q = [0, 2tt] x [0, 2ir]. A 2D random velocity field will be specified 
as the initial condition. The initial fields are initialized by random phase in 
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Fourier spectral space and the initial spectrum is given by [38] 



/ k \ (2s+1) 
E(k, 0) = a s u\k v 1 f — J exp 



i\ (kV 



2 / \ /u^i 



i6) 



where s = 0, 1, 2, • • • , and the normalization constant a s is given by 

a s = (2s + l) s+1 /2 s s!. 

All the results presented below correspond to s = 3, k p = 16, u = {0.1, 0.01} 
and p = 2.7. The lattice size is 512 x 512. The integral length scale L is equal 
to 0.12953. The Reynolds number (Rei = Luo/v) is equal to 111.4. 

In Figs Q-Q, the reconstructed single-particle distribution functions and 
non-equilibrium distribution functions are compared with the real single-particle 
distribution functions and non-equilibrium distribution functions by linear re- 
gression analysis. When uq = 0.1 and t = 1000<5t, it is clear that the recon- 
structed single-particle distribution functions and the non-equilibrium distri- 
bution functions coincide with the real single-particle distribution functions 
and non-equilibrium distribution functions very well for D2Q9 and D2Q17 in 
Figs (JTJ)- (J2J) . The corresponding relative errors E(/ i; fi) are about 0.242% and 
0.194%, respectively. The relative errors E(f < > ncq \ f} neq ^) are about 16.735% 
and 15.782% for the single-particle non-equilibrium distribution functions of 



D2Q9 and D2Q17, respectively. If Eq.(60) by Imamura et al [31] is used to 



calculate the single-particle non-equilibrium distribution functions, the rel- 
ative errors E(/f ncq) , if eq) ) are up to about 21.65% and 18.13% for D2Q9 



and D2Q17, respectively. We also adopted Eqs. (61) in |9J and (62)in [32] 
to do the same calculations. The relative errors Fj(f- ncq \ /^ neq ' ) ) of the single- 
particle non-equilibrium distribution functions can be up to about 80% at 
many lattice nodes. In Fig. [5J the numerical relation between // neq ^ and /^ neq ^ 
for the method in [32]. The mean relative error E(// ne , f} neq ^) is larger than 
43.74% for D2Q9. In the statistical procedure, we ignore the points with very 
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small ft q) and f}™* (ff^JT^ < 10 ~ 3 ) for the method in [32]. Here, we 



f(neq) / ^(neq) ?(neq) 

must point out that when f\ n and f}* 10 ^ are very small, the relative errors 
E(/f neq \ f} neqS> ) of the methods in [9"f3"2"] are very large. In such a circumstance, 



the relative error of the non-equilibrium distribution functions by Eq. uTfh is 



also a bit larger, but it still less than that computed by Eq. (60) [31] and 
much less than that computed by Eqs. (61 )~ (62) of [9] and [32], respectively. 
Similar results can be observed for the case of uq = 0.01 at t — 100005i for 
D2Q9 and D2Q17. For the simplicity of presentation, they are not provided 
here. 



In addition we also found that when the single-particle distribution functions 
and non-equilibrium distribution functions are reconstructed, the results from 
D2Q17 model show a better accuracy than that of D2Q9 model. Meanwhile, 
from the both models, more accurate results can be gained when the Mach 
number is reduced. Such results are very reasonable, and can be understood 
as follows. First, D2Q17 model is more accurate to approach Maxwell dis- 
tribution function in discrete velocity spaces than D2Q9 model. Second, low 
Mach number will lead to a reduction of the truncated errors for approach- 
ing Maxwell distributions and a better recovering Navier-Stokes equation. It is 
proved [37] that D2Q17 model can eliminate the third-order term of statistical 
velocity in recovered Navier-Stokes equation. 

Finally, attention is turned to the comparison of vorticity by the real /j(x, t) 
and the reconstructed /«(x, t) in Figs. where the vorticity contour fig- 

ures are given for Uq — 0.1 and Uq = 0.01, respectively. In order to show the 
quantitative sense of the vorticity reconstruction error, we choose 100 and 
1000 time-series samples for uq = 0.1 and uq = 0.001, respectively. The Ir- 
relative departures of the reconstructed vorticity are 0.02% ±0.0014% (D2Q9, 
M = 0.1), 0.005% ± 0.0003% (D2Q17, u = 0.1), 0.01% ± 0.0026% (D2Q9, 
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M = O.Ol) and 0.003% ± 0.0005% (D2Q17, u = 0.01). The agreement is very 
good. 

In all, the proposed two operators can reconstruct the single-particle distri- 
bution functions and non-equilibrium distribution functions accurately and 
effectively. It can be shown that the two reconstruction operators are very 
flexible to apply to other discrete velocity models of lattice Boltzmann equa- 
tion. 



4-2 The rates of convergence 



In order to validate the approach behaviors versus different grid sizes, we give 
the convergence properties of the D2Q9 and D2Q17 models by different mesh 
scales. The 2D Taylor-Green vortex problem is chosen as the intial fields 



u 



Acos ( &i x ) sin ( fc 2 y ) (£) 



A^sva{k 1 x)cos{k 2 y)F{t) 



57) 



V 



Po 



4 



k 2 

cos(2kix) + — |cos(2fc 2 y) 
fc 2 



F 2 (t) 



where F(t) = exp [—v{k\ + k 2 )t\, A = 0.1, k\ = k 2 = 4 and p = p c 2 . The 
computational domain Q = [0,27r] 2 and Re = 10000. The periodic boundary 
conditions are applied in both directions. The initial distribution functions 
are initialized by the reconstruction operator. The reconstruction L 1 and L 2 
relative errors of the distribution functions are calculated at the time steps 
n = {2000, 4000, 6000, 8000, 10000} corresponding to the mesh resolutions h = 
{1/32, 1/64, 1/96, 1/128, 1/160} respectively. In Fig. [8j the relative errors are 
given in the log-log coordinates. From the results, it is clear that for the 
D2Q9 model and the D2Q17 model, they nearly have the same convergence 
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rates which are approximately equal to 2.6. However, the relative errors of 
the D2Q17 model are smaller than that of the D2Q9 model. That means the 
reconstruction precision can be improved when the number of the discrete 



velocity increases. This conclusion is consistent with the result in Sec. 4.1 



4-3 Coupling computations of FVM and LBM for lid-driven cavity flows 

In order to illustrate the feasibility of the recommended reconstruction opera- 
tor, the lid-driven cavity flow is simulated by the coupled LBM-FVM method. 
The computational domain is decomposed in two regions in which the LBM 
and FVM methods are used respectively (see Fig. [9} (a)). The coarseness and 
fineness of the grids can adjusted according to the zone spatial scale in each 
region. If the grid systems at the interface of overlap subregions are not iden- 
tical, space interpolation at the interface is required when transferring the 
information at the interface. In this paper, the identical mesh structures are 
used for FVM and LBM for convenience to avoid the spatial interpolation 
(see Fig. [9]-(b)). In order to implement the coupling computations, the overlap 
Schwartz alternative procedure is used to handle the computations. 

Numerical simulations were carried out for cavity flow of Re = 100, 400 and 
1000 on a grid 200 x 200. The characteristic length of square cavity is L — 1. 
The boundaries of the cavity are stationary walls, except the top-boundary 
with a uniform tangential velocity (u tt R e =wo = 3.33 x 10 -3 , Ut,R e =4oo = 1-33 x 



10 3 , Mt,i? e =iooo = 3.33 x 10 2 ). Fig. 10 shows plots of the stream function 



for the Reynolds number considered. These plots give a clear picture of the 
overall flow pattern and the effect of Reynolds number on the structure of 
the recirculating eddies in the cavity. The smoothness of the stream function 
distribution, especially around the overlap region confirms the correctness of 
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the information transfer at the interface. To further quantify these results, 
the velocity profiles along the vertical and horizontal centerlines of the cavity 



are shown in Fig. [TTJ The results are in close agreement with the benchmark 
solution [39J. The smoothness and consistency of velocity distribution in the 



overlap region is presented in Fig. 12 where a local, enlarged view of the 
vector plot in the overlap region is shown. Clearly, the vectors in the overlap 
region are quiet consistent between the LBM results and the FVM results. 



Figs. 13 and 14 show the contours of horizontal and vertical velocity. It is seen 
that these physical quantities are all smooth across the interface. According 
to the authors' numerical experience, the smoothness of vorticity contour is 
the most difficult to obtain for such coupled computation, because vorticity if 
the derivative of velocity. The contours of vorticity distribution are shown in 



Fig. 15 Over all, the smoothness on the overlap region are quite good, with a 
minor bumpiness of the left-hand vortex contours for the case of Re = 100. 



In all, by the proposed lifting relation, we can couple the mesoscopic LBM 
with FVM to implement the domain decomposition coupling-computations. 
This paves the way for implementing multiscale computations based on LBM 
and macro-numerical methods of finite-family. 



It should be noted that we also tried the coupling computations based on the 



distribution function fi(x,t) reconstructed by Eq. (61) of |Q] and (62) of [32] . 
Unfortunately, all of our tries were unsuccessful and converged solutions could 
not be obtained. 
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Conclusion 

In this paper, we derive the relation to lift the macroscopic variables to the 
microscopic variables for LBM. Two methods of derivation are conducted and 
they lead to the same result. Numerical tests demonstrate that the derived 
lifting relation possesses good precision. The proposed lifting relation offers a 
way to implement the multiscale-computations involving LBM more efficiently 
and robustly. 
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Fitting line slope = 0.83265 
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(b)Linear regression between fj ncq \x,t) and fi Ueq> (x, t) 

Fig. 1. Linear regression (D2Q9, u = 0.1, t = 100(Wi, i = 2): (a)Fit the 
line /i(x,t) = afi(x,t) + 6,where a = 0.99758 and 6 = 0.00135;(b)Fit the line 
/J neq) (x,t) = afl neq) (x,t) + 6,where a = 0.83265 and b = -2.95012 x 10~ 6 . Stan- 
dard deviation: (a)<7 = 0.00308;(b)cr = 9.21597 x 10" 4 . 
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Fitting line slope = 0.99806 
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(b)Linear regression between fj ncq \x,t) and /j inocu (x, £) 

Fig. 2. Linear regression (D2Q17, u = 0.1 t = 1000«5t, i = 2): (a)Fit the 
line fi(x,t) = afi(x,t) + 6,where a = 0.99806 and b = 0.00227;(b)Fit the line 
^ neq) (x,t) = a/f eq) (x,t) + 6,where a = 0.84218 and b = -4.84408 x lO" 6 . Stan- 
dard deviation: (a)<7 = 0.00288;(b)cr = 8.39673 x 10" 4 . 



7(neq) , 



35 




0.57 0.58 0.59 0.60 0.61 



Real fj 

(a)(a)Linear regression between fi(x,t) and fi(x,t) 



0.0020 




-0.002 -0.001 0.000 0.001 0.002 



Real fj neq 

(b)Linear regression between f- ncq \x,t) 

Fig. 3. Linear regression (D2Q9, u = 0.01, t = 10000<W, % = 2): (a)Fit the line 
fi{x,t) = afi(x,t) + 6,where a = 0.99925 and b = 4.41542 x 10" 4 ;(b)Fit the line 
fl neq \x,t) = a/f ncq) (x,t) + 6,where a = 0.83655 and b = -1.51056 x 10~ 8 . Standard 
deviation: (a)<r = 3.52548 x 10" 4 ;(b)o- = 1.01264 x 10~ 4 . 
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(a)Linear regression between fi(x,t) and fi(x,t) 
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(b)Linear regression between f { (x, t) and f} n (x,t) 

Fig. 4. Linear regression (D2Q17, u = 0.01, t = 10000<5i, i = 2): (a)Fit the line 
fi(x,t) = afi(x,t) + 6,where a = 0.99963 and b = 2.17 x 10~ 4 ;(b)Fit the line 
,/J neq) (x,t) = o/f neq) (x,t) + 6,where a = 0.84764 and b = -2.17758 x lO" 8 . Standard 
deviation: (a)o- = 3.37821 x 10~ 4 ;(b)cT = 9.47431 x 10~ 5 . 
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Fig. 5. Linear regression (D2Q9, uq = 0.1, t = 10005t, i = 2): Fit the line 
fi{x,t) = afi(x,t) + 6,where a = 0.48088 and b = -0.248003 x 10~ 6 . 
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(a)Vorticity contour plot by the real fi(x,t) 




(b)Vorticity contour plot by the reconstructed fi(x,t) 



Fig. 6. Vorticity contour plots (D2Q9, uo = 0.1, t = 1000<5t): (a)Vorticity contour 
plot by the real fi(x,t) ; (b) Vorticity contour plot by the reconstructed /j(x, i) 
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(a)Vorticity contour plot by the real fi(x,t) 




(b)Vorticity contour plot by the reconstructed /j(x, t) 

Fig. 7. Vorticity contour plots (D2Q9, u = 0.01, t = lOOOOJt): (a)Vorticity contour 
plot by the real fi(x,t) ; (b) Vorticity contour plot by the reconstructed /j(x, i) 
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(b) the D2Q17 model 
Convergence rates of the reconstruction operator for D2Q9 and D2Q17. 
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(a)Interface structure between two regions of FVM and LBM 
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(b) Grid layout for a 2D lid-driven cavity (200 x 200) 



Fig. 9. Geometric structure and mesh partition: (a)Interface structure between two 
regions of FVM and LBM; (b)Grid layout for a 2D lid-driven cavity (200 x 200) 
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(a) Re = 100 (b) Re = 400 (c) Re = 1000 

Fig. 10. Contour plots of streamline for different Reynolds numbers 
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(b) Vertical velocity profiles 



Fig. 11. Comparisons between Ghia's benchmark solutions and coupling solutions 
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Fig. 12. Enlarge vector plots in overlap regions 
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Fig. 13. Contour plots of horizontal velocity for different Reynolds numbers 
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(a) Re = 100 (b) Re = 400 (c) Re = 1000 

Fig. 14. Contour plots of vertical velocity for different Reynolds numbers 




(a) Re = 100 (b) Re = 400 (c) Re = 1000 

Fig. 15. Contour plots of vorticity for different Reynolds numbers 
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